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I. INTRODUCTION 

Material properties are often controlled by the com- 
plex micro-structures that form during non-equilibrium 
processing. In general terms the dynamics that occur 
during the processing are controlled by the nature and 
interaction of the topological defects that delineate the 



spatial patterns. For example, in spinodal decomposition 
the topological defects are the surfaces that separate re- 
gions of different concentration. The motion of these sur- 
faces is mainly controlled by the local surface curvature 
and non-local interaction with other surfaces or bound- 
aries through the diffusion of concentration. In contrast, 
block-copolymer systems form lamellar or striped phases 
and the topological defects are dislocations in the striped 
lattice. In this instance the defects interact through long- 
range elastic fields. 

One method for modeling the topological defects is 
through the use of 'phase fields'. These can be thought of 
as physically relevant fields (such as concentration, den- 
sity, magnetization, etc.), or simply as auxiliary fields 
constructed to produce the correct topological defect mo- 
tion. In constructing phenomenological models it is often 
convenient to take the former point of view, since physical 
insight or empirical knowledge can be used to construct 
an appropriate mathematical description. In this paper 
the construction of a phase field model 0] to describe 
the dynamics of crystal growth that includes elastic and 
plastic deformations is described. The model differs from 
other phase field approaches J3, 13, iJ] in that the model 
is constructed to produce phase fields that are periodic. 
This is done by introducing a free energy that is a func- 
tional of the local time averaged density field, p(r, t). In 
this description the liquid state is represented by a uni- 
form p and the crystal state is described a p that has the 
same periodic crystal symmetry as a given crystalline lat- 
tice. This description of a crystal has been used in other 
contexts!^ 13, but not previously for describing material 
processing phenomena. For simplicity this model will be 
referred to as the Phase Field Crystal model, or PFC 
model for short. 

This approach exploits the fact that many properties 
of crystals are controlled by elasticity and symmetry. As 
will be discussed in latter sections, any free energy func- 
tional that is minimized by a periodic field naturally in- 
cludes elastic energy and symmetry properties of the pe- 
riodic field. Thus any property of a crystal that is de- 
termined by symmetry (e.g., relationship between elas- 
tic constants, number and type of dislocations, low-angle 



2 



grain boundary energy, coincident site lattices, etc.) is 
also automatically incorporated in the PFC model. 

The purpose of this paper is to introduce and motivate 
this new modeling technique, discuss the basic properties 
of the model and to present several applications to tech- 
nologically important non-equilibrium phenomena. In re- 
mainder of this section a brief introduction to phase field 
modeling techniques for uniform and periodic field is dis- 
cussed and related to the study of generic liquid/crystal 
transitions. In the following section a simplified PFC 
model is presented and the basic properties of the model 
are calculated analytically. This includes calculation of 
the phase diagram, linear elastic constants and the va- 
cancy diffusion constant. 

In Section (jIII|l the PFC model is applied to a number 
of interesting phenomena including, the determination 
of grain boundary energies, liquid phase epitaxial growth 
and material hardness. In each of these cases the phe- 
nomena are studied in some detail and the results are 
compared with standard theoretical results. At the end 
of this section sample simulations of grain growth, crack 
propagation and a reconstructive phase transition are 
presented to illustrate the versatility of the PFC model. 
Finally a summary of the results are presented in Section 



A. Uniform Fields and Elasticity 

Many non-equilibrium phenomena that lead to dy- 
namic spatial patterns can be described by fields that 
are relatively uniform in space, except near interfaces 
where a rapid change in the field occurs. Classic exam- 
ples include order/disorder transitions (where the field 
is the sublattice concentration), spinodal decomposition 
(where the field is concentration) , dendritic growthj^ 
and eutecticsQ. To a large extent the dynamics of these 
phenomena are controlled by the motion and interaction 
of the interfaces. A great deal of work has gone into con- 
structing and solving models that describe both the in- 
terfaces ('sharp interface models') and the fields ('phase 
field models'). Phase field models are constructed by 
considering symmetries and conservation laws and lead 
to a relatively small (or generic) set of sharp interface 
equations 10]. 

To make matters concrete consider the case of spinodal 
decomposition in AlZn. If a high temperature homoge- 
neous mixture of Al and Zn atoms is quenched below the 
critical temperature small Al and Zn rich zones will form 
and coarsen in time. The order parameter field that de- 
scribes this phase transition is the concentration field. To 
describe the phase transition a free energy is postulated 
(i.e., made up) by consideration of symmetries. For spin- 
odal decomposition the free energy is typically written as 
follows, 



dV 



(/(,/.)+i^|WlV2), (1) 



where /(</>) is the bulk free energy and must contain two 
wells one for each phase (i.e., one for Al-rich zones and 
one for Zn-rich zones). The second term in Eq. takes 
into account the fact that gradients in the concentration 
field are energetically unfavorable. This is the term that 
leads to a surface tension (or energy/length) of domain 
walls that separate Al and Zn rich zones. The dynamics 
are postulated to be dissipative and act such that an 
arbitrary initial condition evolves to a lower energy state. 
These general ideas lead to the well known equation of 
motion; 
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where, F is a phenomenological constant. The Gaus- 
sian random variable, rjc, is chosen to recover the correct 
equilibrium fluctuation spectrum, has zero mean and two 
point correlation, 

(r7,(f,i)77,(r',<') = rkBT{V^r5{r- r')5(i - t'). (3) 

The variable a is equal to one if is a conserved field, 
such as concentration, and is equal to zero if is a non- 
conserved field, such as sublattice concentration. 

A great deal of physics is contained in Eq. jSJ and 
many papers have been devoted to the study of this equa- 
tion. While the reader is refereed to for details, the 
only salient points that will be made here is that 1) the 
gradient term and double well structure of /{(/)) in Eq. 

lead to a surface separating different phases and 2) 
the equation of motion of these interfaces is relatively 
independent of the form of f{(f>). For example it is well 
known ^3 that if (j) is non-conserved the normal velocity, 
Vn, of the interface is given by; 



Vn 



A 



(4) 



where k is the local curvature of the interface and A is di- 
rectly proportional to the free energy difference between 
the two phases. If cj) is conserved then the motion of the 
interface is described by the following set of equations 
M, 



At(0) = doK + l3Vn 
dn/dt = JDVV, 



(5) 



where /i = SF/Scj) is the chemical potential, do is the cap- 
illary length, f3 is the kinetic under-cooling coefficient, n 
is a unit vector perpendicular to the interface position, D 
is the bulk diffusion constant and 0^ and 0~ are positions 
just ahead and behind the interface respectively. 

It turns out that Eqs. (@J and ^ always emerge when 
the bulk free energy contains two wells and the local free 
energy increases when gradients in the order parameter 
field are present 0. In this sense Eqs. Q and (0) can 
be thought of as generic or universal features of systems 
that contain domain walls or surfaces. As will be dis- 
cussed in the next subsection a different set of generic 
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features arise when the field prefers to be periodic in 
space. Some generic features of periodic systems are that 
they naturally contain an elastic energy, are anisotropic 
and have defects that are topologically identical to those 
found in crystals. A number of research groups have built 
these 'periodic features' into phase field models describ- 
ing uniform fields. This approach has some appealing 
features, as one can consider mesoscopic length and time 
scales. But it can involve complicated continuum mod- 
els. For example, in Refs. 0,01) a continuum phase-field 
model was constructed to treat the motion of defects, 
as well as their interaction with moving free surfaces. 
Although such an approach gives explicit access to the 
stresses and strains, including the Burger's vector via 
a ghost field, the interactions between the nonuniform 
stresses and plasticity are complicated, since the former 
constitutes a free-boundary problem, while the latter in- 
volves singular contributions to the strain, within the 
continuum formulation. 



form for (f) the free energy becomes, 



- = kA' 

a 

KA 



2 2KA^ 



a" 



(Aa)^ 



dVf (0) 
dVf{4>), (7) 



where Aa = a — ao- At this level of simplification it can 
be seen that the free energy per unit length is minimized 
when a = Got ov ao is the equilibrium periodicity of the 
system. Perhaps more importantly it highlights the fact 
the energy can be written in a Hooke's law form (i.e., 
E = Eo + (fcAa)2) that is so common in elastic phenom- 
ena. Thus a generic feature of periodic systems is that 
for small perturbation (e.g., compressions or expansions) 
away from the equilibrium they behave elastically. This 
feature of periodic systems will be exploited to develop 
models for crystal systems in the next section. 



B. Periodic Systems 

In many physical systems periodic structures 
emerge. Classic examples include block-copolymers 
[Tl| . Abrikosov vortex lattices in superconductors 
and oil/water systems containing surfactants |0| and 
magnetic thin films. In addition many convective 
instabilities such as Rayleigh-Benard convection 

and a Margonoli instability lead to periodic structures 
(although it is not always possible to describe such sys- 
tems using a free energy functional). To construct a free 
energy functional for periodic systems it is important 
to make the somewhat trivial observation that unlike 
uniform systems, these systems are minimized by spatial 
structures that contain spatial gradients. This simple 
observation implies that in a lowest order gradient 
expansion the coefficient of jV^p in the free energy (see 
Eq. (^)) is negative. By itself this term would lead to 
infinite gradients in so that the next order term in the 
gradient expansion must to be included (i.e., jV^t/jp). 
In addition to these two terms a bulk free energy with 
two wells is also needed, so that a generic free energy 
functional that gives rise to periodic structures can be 
written. 



J" 
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where K and Oq are phenomenological constants. 

Insight into the influence of the gradient energy terms 
can be obtained by considering a solution for (p of the 
form (p = ^ sin(27ra;/a). For this particular functional 



C. Liquid/Solid systems 

In a liquid/solid transition the obvious field of interest 
is the density field since it is significantly different in the 
liquid and solid phases. More precisely the density is 
relatively homogeneous in the liquid phase and spatially 
periodic (i.e., crystalline) in the solid phase. The free 
energy functional can then be approximated as; 



dr[H{6p)] 



dr 



(8) 



where / and G are to be determined and 5p is the devi- 
ation of the density from the average density (p) . Under 
constant volume conditions Sp is a conserved field, so that 
the dynamics are given by; 



dSp 
'dt 



, dT 
dSp 



(9) 



where 77 is a Gaussian random variable with zero mean 
and two point correlation. 



(r;(f,<)77(r', i')) = ThV^Sif- r')(5(f - t'). 



(10) 



To determine the precise functional form of the oper- 
ator G(V2) it is useful to consider a simple liquid since 
Sp is small and / can be expanded to lowest order in Sp, 
i.e., 



liq 



fio) + fii)^Sp) + ^{Sp)' 



(11) 



where = {d'f /dSp')sp=o- In this limit Eq. ® takes 
the form 



dSp 
'dt 



-FV' 



/(2)+G(V2) 



(12) 
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FIG. 1: The points correspond to an experimental liquid 
structure factor for ^®Ar at 85°K taken from J. Yarnell, M. 
Katz, R. Wenzel and S. Koenig, Phys. Rev. A, 7, 2130 (1973) 
[ist . The line corresponds to a best fit to Eq. 1221 . 



which can be easily solved to give, 



A typical liquid state structure factor and the corre- 
sponding LUq is shown in Fig. (Q. Thus GCV^) can be 
obtained for any pure material through Eq. H17(l . 

In the solid state the density is unstable to the forma- 
tion of a periodic structure (i.e., to forming a crystalline 
solid phase) and thus ujq must go negative for certain 
values of q. This instability is taken into account by the 
temperature dependence of Z^^-*, i.e.. 



/ 



(2) 



a{T — T„ 



(18) 



Thus when T > Tm, Wq is positive and the density is 
uniform. When T < Tm, Wq is negative (for some values 
of q) and the density is unstable to the formation of a pe- 
riodic structure. To properly describe this state, higher 
order terms in Sp must be included in the expansion of 
f{Sp), since 6p is no longer small. Before discussing the 
properties of a specific choice for f{5p) it is worth point- 
ing out some generic elastic features of such a model. 

As illustrated in the Section Ijl B|) a free energy that is 
minimized by a periodic structure has 'elastic' properties. 
The elastic constants of the system can be obtained by 
formally expanding around an equilibrium state in the 
strain tensor. If the equilibrium state is defined to be 
Spo{r) and the displacement field is u, then dp can be 
written Sp{r) = dpo{f+ u) + e, where e will always be 
chosen to minimize the free energy. Expanding to lowest 
order in the strain tensor gives. 



Sp{q,t) = e-''-'"<5p(g,0) 



r]{lt'), (13) 



where, q is the wavevector, Qq = /'^^ + G{q'^), 6p is the 
Fourier transform of dp, i.e.. 



6p{q,t)= / dfe'?-^6p{f,t)/{2TTf. 



(14) 



and d is the dimension of space. The structure factor, 
S{q,t) = is then; 



S{q,t) = e-^'>'''^^'Siq,0) + ^(l 



-q'QqTt 



(15) 



In a liquid system the density is stable with respect to 
fiuctuations which implies that > 0. The equilibrium 
liquid state structure factor, Sf^^{q) = S{q, oo), then be- 
comes. 



(16) 



This simple calculation indicates the phase field 
method can model a liquid state if the function diq is 



replaced with kBT/Sf^ (q), or 



G{q)^kBT/S2-f<-'l 



(17) 



df a 



where Cij^ki is the elastic constant given by. 



a. 



1 d'^H 



ij,kl 



2! duijUki 



(19) 



(20) 



In Eq. H2U|) the Einstein summation convection is used 
and Uij represents the usual components of the strain 
tensor, i.e.. 



dui 



duj ^ dui dui 



duj dui dui duj 



(21) 



and the subscript o indicates that the derivative are eval- 
uated at 6p = 6po (i.e., Spo = Uij = 0). While Eq. (^01) is 
somewhat formal and difficult to use for a specific model 
it does highlight several important features. Eq. H2(J|I 
shows that the elastic constants are simply related to the 
curvature of the free energy along given strain directions 
Perhaps more importantly Eq. (|2()|l shows that Cij.ki is 
proportional to H which is a function of the equilibrium 
density field Spo- Thus if the free energy is written such 
that is minimized by a Spo that is cubic, tetragonal, 
hexagonal, etc., then Cij^ki will automatically contain 
the symmetry requirements of that particular system. 
In other words the elastic constants will always satisfy 
any symmetry requirement for a particular crystal sym- 
metry since Cij^ki is directly proportional to a function 
that has the correct symmetry. This also applies to the 
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type or kind of defects or dislocations that can occur in 
any particular crystal system, since such deformation are 
determined by symmetry alone. 

In the next section a very simple model of a liq- 
uid/crystal transition will be presented and discussed in 
some detail. This model is constructed by providing the 
simplest possible approximation for f{Sp) that will lead 
to towards a transition from a uniform density state (i.e., 
a liquid) to a periodic density state (i.e., a crystal). 



II. SIMPLE PFC MODEL: BASIC PROPERTIES 

In this section perhaps the simplest possible periodic 
model of a liquid/crystal transition will be presented. 
Several basic features of this model will be approximated 
analytically in the next few subsections. This includes 
calculations of the phase diagram, the elastic constants 
and the vacancy diffusion constant. 

A. Model 

In the preceding section it was shown that a partic- 
ular material can be modeled by incorporating the two 
point correlation function into the free energy through 
Eq. ifTTji . It was also argued that the basic physical fea- 
tures of elasticity are naturally incorporated by any free 
energy that is minimized by a spatially periodic function. 
In this section the simplest possible free energy that pro- 
duces periodic structures will be examined in detail. This 
free energy can be constructed by fitting the following 
functional form for G, 



2\2 



G(V^) = X{qi + V^) 



(22) 



to the first order peak in an experimental structure fac- 
tor. As an example such a fit is shown for argon in Fig. 
^ . At this level of simplification the minimal free energy 
functional is given by; 



aAT + X 



Wo 



u^]. (23) 



In principle other non- linear terms (such as p"^) can be 
included in the expansion but retaining only simplifies 
calculations. The dynamics of p are then described by the 
following equation. 



dp 
dt 



-rv> + r;. = -rv^ 



dp 



(24) 



For convenience it useful to rewrite the free energy in 
dimensionless units, i.e.. 



X = rqo, 



In dimensionless units the free energy becomes 



TXqlt. (25) 



F=^= I dx 

■J~ O 



(26) 



where To = X?q^ ju and 

c^(V2) =r + (l + V2)2. (27) 
The dimensionless equation of motion becomes, 

^=V2(c.(V2)^ + V'^)+C, (28) 

where, (C(^i, ti)C(n3, T2)) = V\I'^6(tx - ri)b{Tx - T2) and 
V = ukBTqi-^l\^. 

Equations H26|l . H27|l and 128|) describe a material with 
specific elastic properties. In the next few sections the 
properties of this 'material' will be discussed in detail. As 
will be shown, some of the properties can be adjusted to 
match a given experimental system and others cannot be 
matched without changing the functional form of the free 
energy. For example the periodicity (or lattice constant) 
can be adjusted since all lengths have been scaled with 
q^. The bulk modulus can be also be easily adjusted 
since the free energy has been scaled with A, u and qo- 
On the other hand this free energy will always produce 
a triangular lattice in two dimensions; 5, 6]- To obtain a 
square lattice a different choice of the non-linear terms 
must be made. This is the most difficult feature to vary as 
there are no systematic methods (known to the authors) 
for determining which functional form will produce which 
crystal symmetry. Cubic symmetry can be obtained by 
replacing i]}^ with |Vi/;|'^ |16l[l7| . 

In the next few subsections the properties of this free 
energy and some minor extensions will be considered in 
one and two dimensions. The three dimensional case will 
be discussed in a future paper. 



B. One dimension 

In one dimension the free energy described by Eq. 1)26(1 
is minimized by a periodic function for small values of 
i/'o and a constant for large values. To determine the 
properties of the periodic state it is useful to make a one 
mode approximation, i.e., "0 « Asm{iqx) -I- -00, which is 
valid in the small r limit. Substitution of this function 
into Eq. gives. 



- = ^ / dx 



_q 

27r 
2 







02 



2 



4 



3^2 



(29) 



where ujq is the fourier transform of u!{W^), i.e., Cuq = 
r + {1 — g2)2_ Minimizing Eq. with respect to q 

gives the selected wavevector, q* — 1. Minimizing !F 
with respect to A gives, A'^ = -4(0;,. /3 -I- i/'o)- This 
solution is only meaningful if A is real, since the density 
is a real field. This implies that periodic solutions only 

r. The minimum free 



exist when r < —S^fZ, since cijq 
energy density is then; 



pP/L = -r2/6 + i/;2(l - r)/2 - B^j^/A. 



(30) 
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FIG. 2: One dimensional phase diagram in the one-mode ap- 
proximation. The solid line is the boundary separating con- 
stant (i.e., liquid) and periodic (i.e., crystal) phases. The 
hatched section of the plot corresponds to regions of liq- 
uid/crystal coexistence. 



Equation (|30|) represents the free energy density of a 
periodic sohition in the one mode approximation. To 
determine the phase diagram this energy must be com- 
pared to that for a constant state (i.e., the state for which 
■0'^ = i/'o) which is, 



(31) 



To obtain the equilibrium states the Maxwell equal 
area construction rule must be satisfied, i.e.. 



4'2 



dij} {ipo) - Aio] = 0, 



(32) 



where V'l is a solution of /ip = /io, V'2 is a solution of 
= Mo and /x(V') = Hp (= ^c) if < {^p > ^c) 
and fi = dJ-/dtpo- Using these conditions it is straight- 
forward to show that for r > —1/4 a periodic state is 
selected for \ipo\ < r/3 and a constant state is se- 
lected when \ipo\ > \f~-rji. For r < —1/4, there can 
exist a coexistence of periodic and constant states. If the 
constant and periodic states are considered to be a liq- 
uid and crystal respectively then this simple free energy 
allows for the coexistence of a liquid and crystal, which 
implies a free surface. The entire phase diagram is shown 
in Fig. (01. 

It is also relatively easy to calculate the elastic energy 
in the one mode approximation. If a = 27r/q is defined 
as the one-dimensional lattice parameter, then the T can 
be written, 



(33) 



where u = (a — ao)/ao is the strain and K is the bulk 
modulus and is equal to; 



K 



(34) 



or for the particular dispersion relationship used here, 
K = — 8(r -I- 3'0o)/3. The existence of such a Hooke's 
law relationships is automatic when a periodic state is 
selected since F always increases when the wavelength 
deviates from the equilibrium wavelength. 



C. Two dimensions 

1. Phase Diagram 

In two dimensions F is minimized by three distinct so- 
lutions for ip. These solutions are periodic in either zero 
dimensions (i.e., a constant), one dimension (i.e., stripes) 
or two dimensions (i.e., triangular distributions of drops 
or 'particles'). The free energy density for the constant 
and stripe solutions are identical to the periodic and con- 
stant solution discussed in the preceding section. The two 
dimensional solution can be written in the general form. 



71, rn 



- V'o, 



(35) 



where, G = nbi + mb2 and the vectors bi and 62 are 
reciprocal lattice vectors. For a triangular lattice the 
reciprocal lattice vectors can be written, 

2n 



bi 



bi = 



aV3/2 
27r 



i/2x + y/2 



where a is the distance between nearest neighbor local 
maxima of ip (which corresponds to the atomic posi- 
tions). In analogy with the one-dimensional calculations 
presented (see Sec. (IIIBIl ) a one mode approximation 
will be made to evaluate the phase diagram and elas- 
tic constants. In a two dimensional triangular system 
a one mode approximation corresponds to retaining all 
fourier components that have the same length. More 
precisely the lowest order harmonics consist of all (n, m) 
pairs such that the vector G has length 271 /{ay/3/2). This 
set of vectors includes {n,m) = (±1,0), (0,±1), (1,-1) 
and (—1,1). Furthermore since ip is a real function the 
fourier coefficients must satisfy the following relation- 
ship, a-a^m = a-ri,m = Ori-m- In addition, by symmetry, 
a±i,o = oo.ii = ai._i = a_i^i. Taking these consider- 
ations into account it is easy to show that in the lowest 
order harmonic expansion for a triangular solution for ijj 
can be represented by; 



i^t = At 



COS (qtx) cos (^qty/Vi 
)S (2qty/Vtj / 



(37) 



7 




r 



FIG. 3; In Fig. (a) the minimum of the free energy is plotted 
as a function of r for tpo = y/—r/2. The solid line is Eq. (I38II 
and the points are from numerical simulations. In Fig. (b) the 
bulk modulus is plotted as a function of r for tpo = \/— r/2. 
The solid line is an analytic calculation {{qtAt)^) and the 
points are from rmmerical simulations. 



Constant 
Phase 




FIG. 4: Two dimensional phase diagram as calculated in one 
mode approximation. Hatched areas in the figure correspond 
to coexistence regions. The small region enclosed by a dashed 
box is superimposed on the argon phase diagram in Fig. 
In this manner the parameter of the free energy functional 
can be chosen to reproduce certain the relevant aspect of a 
liquid /crystal phase transition. 



where At is an unknown constant and qt = 27r/a. Sub- 
stituting Eq. ijSZI) into Eq. minimizing with respect 
to At and qt gives, 

10 V 50^7 2 \ 25 J 

where, 

At = ^ (^7/'o + iV-15r-36V'2j (39) 

qt = ■\/3/2, and 5 is a unit area. The accuracy of this 
one mode approximation was tested by comparison with 
a direct numerical calculation for a range of r's, using 
'Method r as described in Appendix ljVI|) . The time step 
(At) and grid size (Ax) were 0.0075 and 7r/4 respectively 
and a periodic grid of a maximum size of 512Ax x 512Ax 
|18| was used. A comparison of the analytic and numeri- 
cal solutions are shown in Fig. Q for a variety of values 
of r (tpo was set to be y/—r/2). The approximate solution 
is quite close to the numerical one and becomes exact in 
the limit r — > 0. The analytic results can in principle be 
systematically improved by including more harmonics in 
the expansion. 



To determine the phase diagram in two dimensions the 
free energy of the triangular state (i.e., Eq. must be 
compared with the free energy of a striped state (i.e., 
Eq. (|^ ) and a constant state (i.e., Eq. ^^). In ad- 
dition since V' is a conserved field Maxwell's equal area 
construction must be used to determine the coexistence 
regions. The phase diagram arising from these calcula- 
tions is shown in Fig. Q). While this figure does not look 
like a typical liquid/solid phase diagram in the density- 
temperature plane it can be superimposed onto a por- 
tion of an experimental phase diagram. As an example 
the PEG phase diagram is superimposed onto the argon 
phase diagram in Eig. 



2. Elastic Energy 

The elastic properties of the two dimensional trian- 
gular state can be obtained by considering the energy 
costs for deforming the equilibrium state. The free en- 
ergy density associated with bulk, shear and deviatoric 
deformations can be calculated by considering modified 
forms of Eq. i.e., Mx/{1 + C),?//(l + 0) (bulk), 

^/jt{x + (y,y) (shear) and ^/'t(x(l-|-C), y(l — C) (deviatoric). 
In such calculations C represents the dimensionless defor- 
mation, qt = and At is obtained by minimizing F. 
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FIG. 5: The phase diagram of argon. The hatched regions 
correspond to the coexistence regions. The points are from 
the PFC model, phase field model. 

The results of these calculations are given below; 

Fshear/A = F^,,, + a/8 + • • • 
■^deviatoricl ~ -^rnin ^~ 0^/2 "h ' ' ' 

(40) 

where 

a = 4 (s^po + 36i/'2) ^ /75. (41) 

These results can be used to determine the elastic con- 
stants by noting that for a two dimensional system |5lll9|. 

Fbulk = Fmin + [^^H + C'12] + ' ' ' 

Fshear = -Fmi„ + [C44/2] H 

Fdeviatoric = ^'min + ['^H ~ ^"12] + ' ' ' ■ (42) 

The elastic constants are then 

Cii/3 = C12 = C44 = a/4 (43) 

These results are consistent with the symmetries of a two 
dimensional triangular system, i.e., Cn — C12 + 2C44. In 
two dimensions this implies a bulk modulus of B — a/2, 
a shear modules oi fi = a/4, a Poisson's ratio of cr = 1/3, 
and a two dimensional (i.e., F2 = '^Bii/ {B + fi)) Young's 
modulus of Y2 — 2a/ 3. Numerical simulations were con- 
ducted (using the parameters and numerical technique 
discussed in the previous section) to test the validity of 
these approximations for the bulk modulus. The results, 
shown in Fig. Q, indicate that the approximation is 
quite good in the small r limit. 



These calculations highlight the strengths and limita- 
tions of the simphstic model described by Eq. (|^ . On 
the positive side the model contains all the expected elas- 
tic properties (with the correct symmetries) and the elas- 
tic constants can be approximated analytically within a 
one mode analysis. One the negative side, the model as 
written, can only describe a system where Cn — 3Ci2. 
Thus parameters in the free energy can be chosen to pro- 
duce any Cn, but C12 carmot be varied independently. 
To obtain more flexibility a term could be added to 
the free energy. 



3. Dynamics 

The relatively simple dynamical equation for -0 (i-e-, 
Eq. (EHJ) can describe a large number of physical 
phenomena depending on the initial conditions and 
boundary conditions. To illustrate this versatility it 
is useful to consider the growth a crystalline phase 
from a supercooled liquid, since this phenomena simul- 
taneously involves the motion of liquid/crystal inter- 
faces and grain boundaries separating crystals of dif- 
ferent orientations. Numerical simulations were con- 
ducted using the 'Method F as described in Appendix 
(I VI A(l . The parameters for these simulations were, 
{r,iPo,D,Ax,At) = (-1/4, 0.285, 10-^7r/4, 0.0075) on 
a system of size 512Ax x 512Aa; with periodic bound- 
ary conditions. The initial condition consisted of large 
random gaussian fluctuations (amplitude 0.1) covering 
(10x10) grid points in three locations in the simulation 
cell. As shown in Fig. 10) the initial state evolves into 
three crystallites each with a different orientation and a 
well defined liquid/crystal interface. The excess energy 
of the liquid/crystal interfaces is highlighted in Fig. (jSJi) 
where the local free energy density is plotted. 

As time evolves the crystallites impinge and form grain 
boundaries. As can be seen in Fig. © the nature of 
the grain boundary between grains (1) and (3) is signif- 
icantly different from the boundary between grains (2) 
and grain (1) (or (3)). The reason is that the orientation 
of grains (1) and (3) is quite close but significantly differ- 
ent from (2). The low angle grain boundary consists of 
dislocations separated by large distances, while the high 
angle grain boundary consists of many dislocations piled 
together. A more detailed discussion of the grain bound- 
aries will be given in Section Ijlll A|l . Even this small 
sample simulation illustrates the flexibility and power of 
the PFC technique. This simulation incorporates, the 
heterogeneous nucleation of crystallites, crystallites with 
triangular symmetry and elastic constants, crystallites 
of multi orientations, the motion of liquid/crystal inter- 
faces and the creation and motion of grain boundaries. 
While all these features are incorporated in standard mi- 
croscopic simulations (e.g., molecular dynamics) the time 
scales of these simulations are much longer than could be 
achieved using microscopic models. 

One fundamental time scale in the PFC model is the 
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FIG. 6: Heterogeneous nucleation of a three crystallites in a 
supercooled liquid. The grey scale in a) , b) and c) corresponds 
to the density field, and in d), e) and f) to the smoothed 
local free energy. The configurations are taken at times t = 
300,525 and 3975 for a)+d), b)+e) and c)+f) respectively. 
(Note only a portion of the simulation is shown here) 



FIG. 7: Vacancy diffusion times. In this figure the grey scale 
is proportional to the ij)(f,i) — ^eq. The times shown are a) 
t=0, b) t=50, c) t=100 and d) t=150. 



becomes, 

^ = [(c. + 3 + 2^o0t + 0?)) S,p] , (44) 

where (j^t — ''Pt — tpo (see Eq. (|35|) ). The perturbation, 
S'tp, is then expanded as follows, 



iqt{nx+{n+2m)y /\/3)+iQ-f 



(45) 



diffusion time. To envision mass diffusion in the PFC 
model it convenient to consider a perfect equilibrium 
configuration with one 'particle' missing. At the atomic 
level this would correspond to a vacancy in the lattice. 
Phonon vibrations would occasionally cause neighbor- 
ing atoms to hop into the vacancy and eventually the 
vacancy would diffuse throughout the lattice. In the 
PFC model the time scales associated with lattice vi- 
brations are effectively integrated out and all that is left 
is long time mass diffusion. In this instance the den- 
sity at the missing spot will gradually increase as the 
density at neighboring sites slowly decreases. Numer- 
ical simulations of this process are shown in Fig. (jT)) 
using Method I (see App. IVI A|l ) with the parame- 
ters ir,^o,'D,Ax,At) = (1/4, 1/4, 0,7r/4, 0.0075). To 
highlight diffusion of the vacancy, the difference between 
ip{f,t) and a perfect equilibrium state (V'*) is plotted in 
this Fig 0. 

The diffusion constant in this system can be obtained 
by a simple linear stability analysis, or Bloch-Floquet 
analysis, around an equilibrium state. To begin the anal- 
ysis the equation of motion for ijj is linearized around ^/;*, 
i.e., ip = V'*(r) + H{r,t). To first order in Sijj Eq. 



Substituting Eg. into Eq. gives; 



dt 



^ ^ ^ ^n,m^l,p^i—n — l,j- 



■m—p 



3-m 



(46) 



where, a) = r + (1 — kf and kf ^ = [iqt + QY + olii + 
2j)V3. 

To solve Eq. (|46|) a finite number of modes 
are chosen and the eigenvalues are determined. Us- 
ing the modes corresponding to the reciprocal lat- 
tice vectors in the one mode approximation ((m, n) = 
(±1,0),(0,±1), (1,-1), (-1,1) and the (0,0) mode gives 
four eigenvectors that are always negative and thus irrel- 
evant and three eigenvalues that have the form —DQ-. 
The smallest D arises from 6o,o mode and can be deter- 
mined analytically if only this mode is used (the other 
eigenvalues correspond to D w 3,9). Since this is the 
smallest D it determines the diffusion constant in the 
lattice. The solution is 



D = 31/-^ 



1 + 9^2 



(47) 
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FIG. 8: Vacancy difFusion. In this figure the average of the 
standard deviation in the x and y directions is plotted as a 
function of time. 



The accuracy of Eq. H47|l was tested by numerically mea- 
suring the diffusion constant for the simulations shown in 
Fig. {Tj) . In this calculation the envelope of profile of Sijj 
was fit to a gaussian {Ae^^ ) and the standard de- 
viation (cr) was measured. The diffusion constant can 
be obtained by noting that the solution of a diffusion 
equation (i.e., dC/dt = DV^C) is C cx e"'''/-*^*, i.e., 
cr^ = Dt/2. In Fig. © is plotted as a function of 
time and the slope of this curve gives, D « 1.22. This is 
quite close the value predicted by Eq. H47|l which is 1.25. 



III. SIMPLE PFC MODEL: APPLICATIONS 

In this section several application of the PFC model 
will be considered that highlight the flexibility of the 
model. In Sec. Ijlll A|l the energy of a grain boundary 
separating two grains of different orientation is consid- 
ered. The results are compared with the Read-Shockley 
equation|20j and shown to agree quite well for small ori- 
entational mismatch. This calculation, in part, provides 
evidence that the interaction between dislocations is cor- 
rectly captured by the PFC model, since the grain bound- 
ary energy contains a term that is due to the elastic field 
set up by a line of dislocations. In Sec. Ijlll B|) the tech- 
nologically important process of liquid phase epitaxial 
growth is considered. Numerical simulations are con- 
ducted as a function of mismatch strain and show how the 
model naturally produces the buckling instability and nu- 
cleation of dislocations. In Sec. (jlll C|l the yield strength 
of poly- (nano-) crystalline materials is examined. This 
is a phenomena that requires many of the features of the 
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FIG. 9; Schematic of a grain boundary. 

PFC model (i.e., multi orientations, elastic and plastic 
deformations, grain boundaries) that are difficult to in- 
corporate in standard uniform phase field models. The 
yield strength is examined as a function of grain size and 
the reverse Hall-Petch effect is observed. Finally some 
very preliminary numerical simulations are presented in 
Sec. (IIII D|) to demonstrate the versatility of the tech- 
nique. This section includes simulations of grain growth, 
crack propagation and reconstructive phase transitions. 



A. Grain Boundary Energy 

The free energy density of a boundary between two 
grains that differ in orientation is largely controlled by 
geometry. In a finite size two-dimensional system the 
parameters that control this energy are the orientational 
mismatch, 6 and an offset distance A, as shown in Fig. 
@. For small 9, 9 controls the number of dislocations 
per unit length and A controls the average core energy. 
For an infinite grain boundary A becomes irrelevant, un- 
less the distance between dislocation is an integer number 
of lattice constants (and the integer is relatively small). 
Nevertheless it is straightforward to determine a lower 
bound on the grain boundary energy in the small 9 limit, 
by directly relating the dislocation density to 9 and as- 
suming that the dislocation cores can always find the 
minimum energy location. The later assumption restricts 
the calculation to providing a lower bound on the grain 
boundary energy. 

For small 9, Read and Shockley were able to derive 
an expression for the grain boundary energy, assuming 
the dislocation core energy was a constant independent 
of geometry. In two dimensions the energy/length of the 



11 



jrain boundary is 



— -^core 



^fl-lnf^) 
87rd V \ / 



(48) 



where h is the magnitude of the Burger's vector, a is the 
size of the dislocation core, d is the distance between 
dislocations, Y2 is the two dimensional Young's modulus 
and Ecore is the energy/length of the dislocation core. 
To estimate the minimum core energy it is convenient 
to assume the core energy is proportional to the size of 
the core '55, i.e., Ecore = Ba?^ where B is an unknown 
constant. The total energy/length then becomes 



F 
L 



= Ba^ 



1 - In 



27ra 



(49) 



To obtain a lower bound on F/L the unknown parameter 
B is chosen to minimize F/L, i.e., B is chosen to satisfy, 
d{F/L)/da = 0, which gives Ba^ = b'^Y2/16nd. Thus the 
free energy per unit length is. 



F _ PY2 /3 /27ra 



(50) 



Furthermore, from purely geometrical considerations, the 
distance between dislocations is d — a/tan(f?), where 9 
is the orientational mismatch. Finally in the small angle 
limit (tan(0) w 0) Eq. l(5njl reduces to, 



(51) 



where the dislocation core size b was assumed to be equal 
to the lattice constant a. 

To examine the validity of Eq. H51|l the grain bound- 
ary energy was measured as a function of angle. In 
these simulations numerical Method I (see App. Ij VI A|) ) 
was used with the parameter set (r, ^/'o, 2?, Aa;, At) — 
(-4/15, 1/5, 0,7r/4, 0.01). The initial condition was con- 
structed as follows. On a periodic grid of size Lq^ X Ly, 
a triangular solution (i.e., Eq. l|S3) for tp was con- 
structed in one orientation between < x < Lx/4: and 
3X2^/4 < a; < Lx- In the center of the simulation (i.e., 
Lx/4: < X < 3Lx/4:) a triangular solution of a differ- 
ent orientation was constructed. A small slab of super- 
cooled liquid was placed between the two crystals so as 
not to influence the nature of the grain boundary that 
emerged. The systems were then evolved for a time of 
t = 10, 000, after which the grain boundary energy was 
measured. Small portions of sample configurations are 
shown in Fig. 10 for = 5.8° and 6 = 34.2° (the grain 
boundary energy is symmetric around 30°). As expected 
the Read-Shockley description of a grain boundary is con- 
sistent with the small angle configuration. In contrast the 
large angle grain boundary is much more complicated and 
harder to identify individual dislocations. 

The measured grain boundary energy is compared with 
Eq. in Fig. ifTT)) . As expected Eq. provides an 
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FIG. 10: In this figure the grey scale corresponds to the mag- 
nitude of the field i/' for a grain boundary mismatch of = 5.8° 
and 6 — 34.2° in figures a) and b) respectively. In Fig. a) 
squares have been placed at defect sites. 




FIG. 11: In this figure the grain boundary energy is plotted 
as a function of mismatch orientation. The points correspond 
to numerical simulations of the PFC model and the solid line 
corresponds to Eq. 15111 . 



adequate description for small angles but not for large 
angles. The Read-Shockley equation does fit the mea- 
sured result for all 9 reasonably well if coefficients that 
enter the equation are adjusted, as has been observed in 
experiment 0, . This fit is shown in Fig. ifT^ . 
The situation is obviously more complicated in three 
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FIG. 12: In this figure the grain boundary energy of the PFC 
model is compared with experiments on Tinl2H. Lead|^ and 
Copper ;22||. 

dimension since another degree of freedom exists. This 
degree of freedom can be visuahzed by considering taking 
one of the crystals shown in Fig. lO and rotating it 
out of the page. The extra degree of free can lead to 
interesting phenomena, such as coincident site lattices 
that significantly alter the grain boundary energy. The 
PFC model should provide an excellent tool for studying 
such phenomena since it is purely a geometrical effect 
that is naturally incorporated in the PFC approach. 

B. Liquid Phase Epitcixial Growth 

Liquid phase epitaxial growth is a common industrial 
method used to grow thin films that are coherent with 
a substrate. The properties of such films depend on the 
structural integrity of the film. Unfortunately flat defect- 
free heteroepitaxial films of appreciable thickness are of- 
ten difficult to grow due to morphological instabilities in- 
duced by the anisotropic strain arising from the mismatch 
between film and substrate lattice constants. Conse- 
quently, there has been a tremendous amount of scientific 
effort devoted to understanding the morphological stabil- 
ity of epitaxially grown filmsjlljl, 0, "2^, '2?^, "2^, "2?, 

The stability and resulting structural properties of epi- 
taxial films are often compromised by at least two dis- 
tinct processes which reduce the anisotropic strain. In 
one process, small mounds or ridges form as the sur- 
face buckles or corrugates to reduce the overall strain 
in the film. This instability to buckling can be predicted 
by considering the linear stability of an anisotropically 
strained film as done by Asaro and Tiller and Grin- 



FIG. 13: Epitaxial growth. Figures a), b), c) and d) corre- 
spond to times t=150, 300, 450 and 600 respectively. The 
grey-scale is proportional to the local density (i.e., ip) in the 
film and the liquid. The substrate is highlighted by a darker 
grey background. To highlight nucleated dislocations, small 
white dots were placed on atoms near the two dislocation 
cores that appear in this configuration. 

feld |23| • The initial length scale of the buckling is deter- 
mined by a competition between the reduction in overall 
elastic energy which prefers mounds and surface tension 
and gravity both of which favor a flat interface. Another 
mechanism that reduces strain is the nucleation of misfit 
dislocations which can occur when the energy of a dislo- 
cation loop is comparable with the elastic energy of the 
strained film. Matthews and Blakeslee ^3 ^^'^ many 
others HHH El HHM El have used various argu- 
ments to provide an expression for the critical height at 
which a flat epitaxially grown fllm will nucleate misfit 
dislocations. 

The two mechanisms are often considered separately 
but it is clear that surface buckling can strongly influ- 
ence the nucleation of misfit dislocations. Typically, as 
the film begins to grow, it will deform coherently by 
the Asaro-Tiller-Grinfeld instability. This leads initially 
to a roughly sinusoidal film thickness with a periodic- 
ity close to the most unstable mode in a linear analysis. 
As time increases, the sinusoidal pattern grows in am- 
plitude and develops cusps or local regions of high cur- 
vature El 0, EE EH with a periodicity similar to that 
of the initial instability although some coarsening may 
occur 0, EEEll- Eventually, the stress at the cusps be- 
come too large and a periodic array of misfit dislocations 
appear which reduces the roughness of the film. These 
dislocations eventually climb to the film/substrate inter- 
face. 

The purpose of this section is to illustrate how the 
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PFC method can be exploited to examine surface buck- 
ling and dislocation nucleation in liquid phase epitaxial 
growth. Modeling this process requires a slight modifi- 
cation of the model to incorporate a substrate that has 
a different lattice constant than the growing film. This 
can be accomplished by changing the operator uj given in 
Eq. ^ to be, 

uj^r+{q^ + V^f, (52) 

where the parameter q controls the lattice constant of 
the growing film and is set to one in the substrate. To 
incorporate a constant mass flux the field ■0 was fixed to 
be, ■0^ at a constant distance {L = lOOAx) above the 
film. 

Numerical simulations were conducted using Method 
I (see App. (|VI All ) for the parameters (r, tpi, Ax, At) = 
(-1/4, .282, .785,0.0075). The width of the film grown 
was Lx — 8 192 Ax, corresponding to a width of roughly 
900 particles. The initial condition was such that eight 
layers of substrate atoms resided at the bottom of the 
simulation cell with a supercooled (r — — 1/4,?/;^ — .282) 
liquid above it. A small portion of a simulation is shown 
in Fig. (|13|l . for q = 0.93. As can be seen in this figure, 
and in Fig. H14() the film initially grows in a uniform fash- 
ion before becoming unstable to a buckling or mounding 
instability. The film then nucleates dislocations in the 
valleys where the stress is the largest. After the disloca- 
tions nucleate the liquid/film interface grows in a more 
regular fashion. To highlight the local elastic energy, the 
free energy is plotted in Fig. (|14|l . As can be seen in 
this figure, elastic energy builds up in the valleys during 
the buckling instability and is released when dislocations 
appear. The behavior of the liquid/film interface was 
monitored by calculating the average interface height and 
width. Both quantities are plotted in Fig. H15|) . The data 
shown in this figure is representative of all simulations 
conducted at different mismatch strains, but the precise 
details varied from run to run. In all cases the width ini- 
tially fluctuates around a*/2 (where a* is the thickness of 
a film layer) during the 'step by step' growth. The aver- 
age width then increases during buckling and decreases 
when dislocations nucleate. While these quantities are 
difficult to measure in situ there is experimental evidence 
for this behavior in SiGe/Si heterostructures [43l |. 

Assigning a value to the critical height. He at which 
dislocations nucleate is very subjective. Typically a first 
wave of dislocations is nucleated at a density that is de- 
termined by the buckling instability. Since this is not the 
correct density to reduce the strain to zero a subsequent 
buckling and dislocation occurs above the first wave. To 
complicate manners the nucleated dislocations climb to- 
wards the substrate/film interface. To illustrate these 
points the dynamics of a sample distribution of defects 
is shown as function of height in Fig. (|16|) . As can be 
seen in this figure the first wave of dislocations appears 
roughly between a film height of six and thirteen layers. 
Comparison of Figs. HlGb l and IjKit i') shows that as time 
evolves the overall distribution of dislocation climbs to- 




FIG. 14: Epitaxial growth. Figures a), b), c), d) and e) 
correspond to times t=150, 300, 450, 600 and 750 respectively. 
The grey-scale is proportional to the free energy density. To 
highlight the excess strain energy in the film the grey scale 
near the defect was saturated. The region enclosed by dashed 
lines corresponds to the configuration shown in Fig. 11311 . 
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Lime 

FIG. 15: Epitaxial growth. In Figs, a) and b) the average 
film/liquid interface height and width is shown as a function of 
time. Both the width and height have been scaled by a* — 27r, 
which is the one mode approximation for the distance between 
layers in the appropriate direction 



14 




5 10 15 5 10 15 -0 1 -0,05 0.05 1 

Hieght/a* Hieght/a* ^ 



FIG. 16: Epitaxial growth. In this figure a histogram of the 
number of defects is shown as a function of height above 
the substrate. Figs, a), b), c) and d) correspond to t — 
300, 450, 600 and 1000. 



FIG. 17: Epitaxial growth. In this figure the He is the critical 
height as defined in the text and e is the mismatch strain 
between the film and substrate. 



ward the surface. To obtain an operational definition of 
He, the average height, H{t) of the first wave of disloca- 
tions was monitored as a function of time. Typically H{t) 
is a maximum when all dislocations in the first wave have 
appeared and then decreases as the dislocation climb to 
the substrate/film interface. He was defined as the max- 
imum value of H{t). 

The critical height, as defined in the preceding para- 
graph, was calculated as a function of mismatch strain, 
e = {a film - a-substrate) / dsubstrate- The equilibrium lat- 
tice constant in the film a fum was obtained by assuming 
it was directly proportional to 1/g (note, in the one mode 
approximation, a = 27r/(y^3)q/2)) and determining the 
constant of proportionality by interpolating to where the 
critical height diverges. The numerical data was com- 
pared with the functional form proposed by Matthews 
and Blakeslee j4^ , i.e.. 



TJc oc i ( 1 + logio 



Hr 



(53) 



in Fig. H17() . This comparison indicates that the data 
is consistent with a linear relationship between e and 
[1 + \og{Hc/a*)]/ {Hc/a*), where the constant of propor- 
tionality depends on whether a compressive or tensile 
load is applied to the substrate. 



C. Material Hardness 

It is well known that mechanical properties of materi- 
als depend crucially on the microstructure and grain size. 



For example. Hall and Fetch |4j, |43 calculated that for 
large grain sizes, the yield strength of a material is in- 
versely proportional to the square of the average grain 
radius. This result is due to the pileup of dislocations at 
grain boundaries and has been verified in many materi- 
als including Fe alloyspl El El, Ni Ni-P alloys 
[sof. Cu j5lj and Pd |5lj. However, for very small grain 
sizes the Hall-Petch relationship must break down, since 
the yield strength cannot diverge. Experimentally it is 
found that materials "soften" at very small grain sizes, 
such that the yield strength begins to decrease when the 
grain sizes become of the order of tens of nanometers. 
This 'inverse' Hall-Petch behavior has been observed in 
in Ni-P alloys Cu and Pd ,51j and molecular dy- 
namics experiments [s^, 01 • Determining the precise the 
crossover length scale and mechanisms of material break- 
down has become increasingly important in technological 
processes as interest in nano- crystalline materials (and 
nano-technology in general) increases. 

The purpose of this section is to demonstrate how the 
PFC approach can be used to study the influence of 
grain size on material strength. In these simulations a 
poly-crystalline sample was created by heterogenous nu- 
cleation (see Sec. (|III D ip for details) in a system with 
periodic boundary conditions in both the x and y di- 
rections. A small coexisting liquid boundary of width 
200Aa: was included on either side of the sample. To ap- 
ply a strain the particles near the liquid/crystal boundary 
(i.e., within a distance of 16 Ax) were 'pulled' by coupling 
these particles to a moving field that fixed the particle 
positions. Initially the system was equilibrated for a total 
time of 4000 (2000 before the field was applied and 2000 
after). An increasing strain was modeled by moving the 




FIG. 18: In this figure the grey scale corresponds to the local 
energy density before a strain is applied. The dark black 
regions on the left and right of the figure are the regions that 
are coupled to the external field. 



field every so many time steps in such a manner that the 
size of the polycrystal increased by 2 Ax. To facilitate re- 
laxation, ip was extrapolated to the new size after every 
movement of the external field. The parameters of the 
simulations to follow were (r, ipsoh Ly, Ax, At) = 

(-0.3, 0.312, 2048Ax,2048Ax, 0.377, 0.79, 0.05) and the 
pseudo-spectral numerical method described in App. 
(|VIB|l was used. 

A sample initial configuration is shown in Fig. (|18|l . 
This particular sample contains approximately 100 grains 
with an average grain radius of 35 particles. As can be 
seen in this figure there exists a large variety (i.e., dis- 
tribution of mismatch orientations) of grain boundaries 
as would exist in a realistic poly-crystalline sample. The 
same configuration is shown after it has been stretched 
in the x direction in Figs. H19|l and (|20|l corresponding 
to strains of 1.9% and 7.8% respectively. 

As the poly-crystalline sample is pulled the total free 
energy was monitored and used to calculate the stress, 
i.e., stress = (PF/dC,"^, where C, is the relative change in 
the width of the crystal. Stress-strain curves are shown 
in Fig. (|21|) as a function of grain size and strain rate. 
In all cases the stress is initially a linear of function of 
strain until plastic deformation occurs and the slope of 
the stress-strain curve decreases. In Fig. (|21b ') the in- 
fluence of strain rate is examined for the initial config- 
uration shown in Fig. (|18|l . It is clear from this figure 
that strain rate plays an strong role in determining the 
maximum stress that a sample can reach, or the yield 
stress, as has been observed in experiments j55| . The 
yield strength increases as the strain rate increases as 
would be expected. 
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FIG. 19: This figure is the same as Fig. IjlH^ except at a 
strain of 1.9%. 




FIG. 20: This figure is the same as Fig. 1181 . except at at 
strain of 7.8%. 



The influence of grain size on the stress-strain relation- 
ship is shown in Fig. H21b ) for four grain sizes. The initial 
slope of the stress-strain curve (which will be denoted Yq 
in what follows) increases with increasing grain size as 
does the maximum stress, or yield stress, sustained by 
the sample. The yield strength and elastic moduh (Iq) 
are plotted as a function of inverse grain size in Figs. 
(|22|l and H23|l respectively for several strain rates. For 
each strain rate the yield stress is seen to be inversely 
proportional to the square root of the average grain size. 
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FIG. 21: In Fig. a) the stress is plotted as a function of strain 
for a system with an average grain radius of 35 particles. The 
solid lines from top to bottom in a) correspond to strain rates 
of 24 X 10"^, 12 X 10"® and 6 x 10"® respectively. In Fig. b) 
the solid lines from top to bottom correspond to systems with 
average grain sizes of 70, 50, 35 and 18 particles respectively. 
In both a) and b) the dashed line corresponds to a unit slope. 



FIG. 22: In this figure the yield stress is plotted as function 
of average grain radius. The solid, empty and starred points 
correspond to strain rates of of 24 x 10~®, 12 x 10~® and 
6 X 10~® respectively. The dashed lines are guides to the eye. 



however the constant of proportionahty decreases with 
decreasing strain rate. Thus the PFC approach is able to 
reproduce the inverse Hall-Petch effect or the softening 
of nano-crystalline materials. 

It would be interesting to observe the cross-over to 
the normal Hall-Petch effect where the yield stress de- 
creases with increasing grain size. However, it is im- 
portant to note that the initial conditions in these sim- 
ulations was set up to explicitly remove the Hall-Petch 
mechanism, i.e., each nano-crystal was defect free. In 
addition thermal fluctuations wore not included in the 
simulations. Nevertheless it is unclear whc;ther or not a 
crossover may occur, due to the fact that low angle grain 
boundaries may act as sources of movable dislocations. 
Further study of this interesting phenomena for larger 
grain sizes would be of great interest. 



D. Other Phenomena 




* 



0.2 - * 
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FIG. 23: The clastic moduli Yq (sec text) is plotted as a 
function grain radius in this figure. The solid, empty and 
starred points correspond to strain rates of of 24 x 10"'', 12 x 
10"® and 6 x 10"® respectively. 



There are many phenomena that the PFC method can 
be used to explore. To illustrate this a few small simula- 
tions were conducted to examine a number of interesting 
phenomena of current interest. In the next few sections 
some preliminary results are shown for grain growth, 
crack propagation and reconstructive phase transitions. 



1. Grain Growth 

When a liquid is supercooled just below the melt- 
ing temperature small crystallites can nucleate homoge- 
neously or heterogencously. The crystallites will grow 
and impinge on neighboring crystallites forming grain 
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boundaries. Depending on the temperature and aver- 
age concentration the final state (i.e., in the infinite time 
limit) may be a single crystal or a coexistence of liq- 
uid and crystal phase since there exists a miscibility gap 
in density for some regions of the phase diagrams. For 
deep temperatures quenches the liquid is unstable to the 
formation of a solid phase and initially an amorphous 
sample is created very rapidly which will evolve into a 
poly-crystalline sample and eventually become a single 
crystal (in the infinite time limit). All these phenomena 
can be studied with the simple PFC model considered in 
this paper. 

In this section the PFC model is used to examine the 
heterogenous nucleation of a poly-crystalline sample from 
a supercooled liquid state. A simulation containing fifty 
initial seeds (or nucleation sites) was conducted. The ini- 
tial seeds were identical to those described in Sec. (Ill C 3|) 
as were all other relevant parameters. The results of the 
simulations are shown in Fig. H24|l . Comparison of Figs. 
(f^lb ) and if^lh ') shows that there is a wide distribution of 
grain boundaries each with a different density of disloca- 
tions (which appear as black dots in the figure) . Compar- 
ison of 1)241) with later configurations indicates that the 
low angle grain boundaries disappear much more rapidly 
than the large angle ones. The simple reason is that it is 
easy for one or two dislocations to glide in such a manner 
as to reduce the overall energy (this is usually accompa- 
nied with some grain rotation). The simulation was run 
for up to a time of t = 50,000 (or approximately 1,200 
diffusion times) and contained approximately 15,000 par- 
ticles. The simulation took roughly 70 hours of cpu on a 
single alpha chip processor (xplOOO). 

2. Crack Propagation 

The PFC model can be used to study the propagation 
of a crack in ductile (but not brittle) material. To illus- 
trate this phenomena a preliminary simulation was con- 
ducted on a periodic system of size (4096Ax, 1024Aa;) for 
the parameters (r, V'o, Ax, At) = (-1.0, 0.49, 7r/3, .05). 
Initially a defect free crystal was set up in the simula- 
tion cell that had no strain in the x direction and a 10% 
strain in the y direction. A notch of size 20Ax x lOAx 
was cut out of the center of the simulation cell and re- 
placed with a coexisting liquid {ij} = 0.79). The notch 
provides a nucleating cite for a crack to start propagat- 
ing. A sample simulation is shown in Fig. H25|) . 

3. Reconstructive Phase Transitions 

The simple PFC model can be used to study a phase 
transition from a state with square symmetry to one with 
triangular symmetry. In the model described by Eq. 
(|26|l a state with square symmetry is metastable, i.e., a 
state with square symmetry will remain unchanged unless 
boundary conduction or fluctuations are present. Bound- 




FIG. 24: Heterogenous nucleation and grain growth. In this 
figure the grey scale corresponds to the smoothed local free 
energy. The figures a), b), c), d), e) and f) correspond to 
times 50, 200, 1000, 3000, 15000, 50000 respectively. 



ary conduction or fluctuations allow for the nucleation of 
a lower energy state which in this particular model is the 
state of triangular symmetry discussed in Sec. I|II C 1} . 
A small simulation was performed to illustrate this phe- 
nomena. In this simulation a crystallite with square 
symmetry coexisting with a liquid was created as an ini- 
tial condition. The parameters for this simulation were 
{r,TPuq,^soi,^x,M) = (1.0, 0.68, 0.52, 1.0, .02). The 
simulations depicted in Fig. 1)26(1 . show the spontaneous 
transition from square structure to the triangular struc- 
ture. Two variants of the triangular structure (differing 
by a rotation of 30 degrees) form in the new phase as 
highhghted in Fig. l(Ml). 

A better method for studying this phenomena is to 
create a free energy that contains both square and trian- 
gular symmetry equilibrium states. This can be done by 
including a |VV'|^ term (which favors square symmetry) 
in the free energy. This is, unfortunately not the most 
convenient term for numerical simulations. A better ap- 
proach is to simply couple two fields in the appropriate 
manner as was done in an earlier publication ,1J. In ei- 
ther case an initial poly-crystalline state can be created 
of one crystal symmetry. 
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FIG. 25: In this figure a portion of a simulation is shown 
where the grey scale corresponds to the local energy density. 
The size of both figures is 2048Aa;x 1024Aa;, where Ax = tt/3. 
Figs, a) and b) are at times t — 25, 000 and 65, 000 after the 
rip was initiated respectively. 



B ■• ■■ il- 




FIG. 26: In this figure the grey scale corresponds the field tp- 
Figs, a), b), c) and d) correspond to times, t = 2, 20, 40 and 
180 respectively. In Fig. d) the solid lines are guides to the 
eye. 



phenomena in three dimensions. 
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VI. APPENDIX: NUMERICAL METHODS 

Equation H28() was numerically solved using the two dif- 
ferent methods work as described below. In what follows 
the subscripts n, i and j are integers that corresponds to 
the number of time steps and distance along the x and y 
directions of the the lattice respectively. Time and space 
units are recovered by the simple relations, t = nAt, 
X = iIS.x and y — jAx. 



A. Method I 

In method I an Euler discretization scheme was used 
for the time derivative and the 'spherical laplacian' ap- 
proximation was used to calculate all Laplacians. For 
this method the discrete dynamics read, 



where fJ-n,i,j is the chemical potential given by; 

^in,^,J = (r + (1 + S/^)^)^n,^.J ~ ^P^.^j. 

All Laplacians were evaluated as follows. 



(54) 



(55) 



^ fn,i,j — |^(/n,z+l,J fn.i — l.j ~^ /n.i.j + l 
+ /n,i,j-l)/2 + (/n,i+l,j + l + /n,^-l,j + l 

+ + /„,^-i,j+i)/4 - 3/„,,j)/(Aa;)2.(56) 



Method II 



IV. SUMMARY 

The purpose of this was paper was to introduce the 
PEG method of studying non-equilibrium phenomena in- 
volving elastic and plastic deformations and then to show 
how the technique can be applied to many phenomena. 
Those phenomena included epitaxial growth, material 
hardness, grain growth, reconstructive phase transitions, 
crack propagation, and spinodal decomposition. In the 
future, we intend to extend this model to study these 



In method II an Euler algorithm was again used for 
the time step, except that a simplifying assumption was 
made to evaluate (r -I- (1 -t- V^)^)'0n,i j) in Fourier space. 
In this approach the fourier transform of ipn.ij was nu- 
merically calculated then multiplied by w(q) and then 
an inverse fourier transform was numerically evaluated 
to obtain an approximation to (r + (1 -|- V^)^)^/'„^i j ). If 
w{q) is chosen to be ■w{q) — r + (1 — then, to within 
numerical accuracy, there is no approximation. In this 
work 'w{q) was chosen to be r -I- (1 — q^Y if w{q) < —2.5 
and w{q) = —2.5 otherwise. Thus w{q) is identical to the 
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exact result for wavevcctors close to g = 1, i.e., the wave- modes arise from the largest negative values ofw(q). This 
lengths of interest. The advantage of introducing a large allows the use of much larger time steps. Other than this 
wavevector cutoff is that the most numerically unstable approximation the method is identical to Method I. 
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